18 November 2020

Intro

Agenda

  • Intro
  • One-sample methods:
    • Pre/post means
    • Interrupted time series
  • Two-sample methods:
    • Treatment/comparison means
    • Comparative interrupted time series
    • Difference in differences
  • Regression methods
    • Fixed effects
    • Linear probabilty models

About me, about this talk

  • Not an authority
  • Been thinking about this stuff
  • It’s been weird having a biostats background and then jumping in with all these economist-types.
  • They probably assumed I was 100% on board with all this econometrics stuff when they asked me to do this HAHAHAAHAa

Estimate, estimator, estimand

Notation

\[ \begin{align} Y(t) ~~&|~~ \mbox{Observed outcome at time } t \\ D = 0 ~~&|~~ \mbox{Control} \\ D = 1 ~~&|~~ \mbox{Treated} \\ t = 1, \ldots, T_0 ~~&|~~ \mbox{Pre-intervention times} \\ t = T_0 + 1, \ldots, T ~~&|~~ \mbox{Post-intervention times} \\ Y^d(t) ~~&|~~ \mbox{Potential outcome with treatment } D=d \mbox{ at time } t \end{align} \]

Counterfactual framework

Two parallel universes

  • one where the intervention occurred
  • one where the intervention did not occur



Policy effect = difference between the outcomes in those two universes

Average treatment effect on the treated (ATT)

Estimand: \[ \begin{align} ATT =& E[Y^1 - Y^0 | D = 1] \\ =& E[Y^1 | D = 1] - E[Y^0 | D = 1] \end{align} \]

Estimator:

\[ \begin{align} ATT &\approx E[Y^1 | D = 1] - \mbox{(some surrogate for } E[Y^0 | D = 1]) \end{align} \]

Application: Terraforming LV-426

Hadley’s Hope

LV-426

Building Better Worlds

The Weyland-Yutani workplace safety program

  • Terraforming workers at the Hadley’s Hope colony on exomoon LV-426 sometimes experience serious workplace injuries.
  • These injuries can lead to expensive claims against the company and delays in the colonization schedule.
  • On October 1, 2176, the Weyland-Yutani corporation implemented a safety program that included a mandatory 3-hour training workshop and disciplinary measures for workers who don’t follow safety protocols.

Data

## # A tibble: 6 x 8
##   colony         time     t colony_effect weather_index treat  post     y
##   <chr>         <int> <int>         <dbl>         <dbl> <dbl> <dbl> <dbl>
## 1 Hadley's Hope     1   -30         -10.4          3.48     1     0  516.
## 2 Hadley's Hope     2   -29         -10.4          3.79     1     0  533.
## 3 Hadley's Hope     3   -28         -10.4          3.36     1     0  556.
## 4 Hadley's Hope     4   -27         -10.4          6.25     1     0  548.
## 5 Hadley's Hope     5   -26         -10.4          8.69     1     0  575.
## 6 Hadley's Hope     6   -25         -10.4          7.94     1     0  578.

Broad research question

What effect does this policy intervention have on monthly operating costs?

Some strategies

  • Pre/post test
  • Interrupted time series (ITS)
  • Comparative ITS
  • Difference in differences

Comparing pre/post means

Operating costs over time: Panel data

Compare means

ATT

\[ \begin{align} ATT =& E[Y^1 - Y^0 | D = 1] \\ =& E[Y^1 | D = 1] - E[Y^0 | D = 1] \end{align} \] \(E[Y^0 | D = 1]\) is in the other universe, so assume

\[ \begin{align} E[Y^0 | D = 1] &\approx E[Y^1(t \le T_0) | D = 1] \end{align} \]

Counterfactual

ATT

Then

\[ \begin{align} ATT &= E[Y^1 - Y^0 | D = 1] \\ &= E[Y^1 | D = 1] - E[Y^0 | D = 1] \\ &= E[Y^1(t > T_0) | D = 1] - E[Y^0(t > T_0) | D = 1] \\ &= E[Y^1(t > T_0) | D = 1] - E[Y^1(t < T_0) | D = 1] \\ & = \bar{X}_{t > T_0} - \bar{X}_{t < T_0} \end{align} \]

t-test

\[ t = \frac{\bar{X}_{post} - \bar{X}_{pre}}{s \sqrt{\frac{2}{n}}} \]

mean(hh$y[hh$post==1]) - mean(hh$y[hh$post==0])
with(hh, t.test(y ~ post))
## [1] 31.31477
## 
##  Welch Two Sample t-test
## 
## data:  y by post
## t = -3.7687, df = 57.917, p-value = 0.0003867
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -47.94785 -14.68169
## sample estimates:
## mean in group 0 mean in group 1 
##        583.1383        614.4531

OLS

Similarly, \[ Y = \beta_0 + \beta_1 \mbox{POST} + \mathbf{\beta_x X} + \varepsilon \]

lm_prepost <- lm(y ~ post + weather_index,
           data=hh)
## # A tibble: 3 x 5
##   term          estimate std.error statistic p.value
##   <chr>            <dbl>     <dbl>     <dbl> <chr>  
## 1 (Intercept)     553.        9.42     58.7  <0.001 
## 2 post             38.5       7.71      4.99 <0.001 
## 3 weather_index     5.11      1.34      3.83 <0.001

Assumptions

  • Outcomes would stay the same over time in the absence of the intervention
  • No unmeasured variables affecting the outcome

Interrupted time series

Operating costs over time

Change in level

Change in level

\[ \begin{align} \mbox{Compare means:} ~~|&~~ Y = \beta_0 + \beta_1 POST + \varepsilon \\ \mbox{ITS (level):} ~~|&~~ Y_{t} = \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \varepsilon_t\\ \end{align} \]

\[ \begin{align} ATT &= E[Y^1 - Y^0 | D = 1] \\ &= E[Y^1 | D = 1] - E[Y^0 | D = 1] \\ & = \hat{\beta}_2 \end{align} \]

Linear model

lm_its <- lm(y ~ t + post,
             data = hh)
## # A tibble: 3 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl> <chr>  
## 1 (Intercept)   610.       8.10      75.3  <0.001 
## 2 t               1.80     0.421      4.27 <0.001 
## 3 post          -22.6     14.6       -1.55 0.127

Change in level and slope

Change in level and slope

Linear model

\[ \begin{align} \mbox{Compare means:} ~~|&~~ Y = \beta_0 + \beta_1 POST + \varepsilon \\ \mbox{ITS (level):} ~~|&~~ Y_{t} = \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \varepsilon_t\\ \mbox{ITS (level & slope):} ~~|&~~ Y_{t} = \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \beta_3 \mbox{TIME} \times \mbox{POST} + \varepsilon_t\\ \end{align} \]

ATT

\[ Y_{t} = \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \beta_3 \mbox{TIME} \times \mbox{POST} + \varepsilon_t \]

So for any \(t \in T \ge T_0\),

\[ \begin{align} ATT(t) &= E[Y^1(t) - Y^0(t) | D = 1] \\ &= E[Y^1(t) | D = 1] - E[Y^0(t) | D = 1] \\ & = \hat{\beta}_2 + (\hat{\beta}_3 \times t) \end{align} \]

Change in level and slope

Linear model

lm_its2 <- lm(y ~ t*post,
             data = hh)
## # A tibble: 4 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl> <chr>  
## 1 (Intercept)   620.       9.73      63.7  <0.001 
## 2 t               2.46     0.557      4.41 <0.001 
## 3 post          -20.4     14.4       -1.42 0.162  
## 4 t:post         -1.47     0.831     -1.77 0.081

Interrupted time series

Assumes the observations from before the intervention are a reasonable surrogate for the counterfactual condition of the test group to represent what would have happened without an intervention.

Accounts for outcomes that were already getting better or worse

Does not account for anything unrelated to the intervention that changed outcomes in the post-intervention period

Exogenous variables?

P-5000 Powered Work Loader

What if there’s some external force causing costs to fluctuate over time?

Comparative methods

Fiorina 161

Fiorina 161, Outer Veil

Fiorina 161 and Hadley’s Hope

Two-group post-mean comparison

Two-group post-mean comparison

Two-group post-mean comparison

\[ t = \frac{\bar{X}_{test} - \bar{X}_{comparison}}{s \sqrt{\frac{2}{n}}} \]

mean(hhf$y[hhf$treat==1]) - mean(hhf$y[hhf$treat==0])
## [1] -44.78596
with(hhf, t.test(y ~ post))
## 
##  Welch Two Sample t-test
## 
## data:  y by post
## t = -7.4716, df = 103.55, p-value = 2.617e-11
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -78.28028 -45.44163
## sample estimates:
## mean in group 0 mean in group 1 
##        590.7673        652.6282

Two-group post-mean comparison

\[ \begin{align} \mbox{Test/comparison post-means:} ~~|&~~ Y = \beta_0 + \beta_1 \mbox{TREAT} + \mathbf{\beta_x X} + \varepsilon \end{align} \]

lm_2grp <- lm(y ~ treat + weather_index,
           data=hhf %>% filter(post==1))
## 
## 
## |term          |  estimate| std.error| statistic|p.value |
## |:-------------|---------:|---------:|---------:|:-------|
## |(Intercept)   | 661.96674|  8.504821| 77.834292|<0.001  |
## |treat         | -75.80414|  7.861046| -9.643010|<0.001  |
## |weather_index |   6.30736|  1.408055|  4.479483|<0.001  |

Two-group post-mean comparison

  • Assumes that the comparison group is a good representation of the counterfactual \(E[Y^0 | D=1]\); specifically, that outcomes for the two groups would be the same if neither or both had received the intervention

  • Accounts for anything unrelated to the intervention that changed outcomes equally for both groups

  • Does not account for any baseline differences in levels or slopes between the two groups or non-intervention shocks that affected only one group

Comparative interrupted time series

Intuition

\[ \begin{align} \mbox{ITS (level):} ~~|&~~ Y_{t} = \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \varepsilon_t\\ \mbox{ITS (level & slope):} ~~|&~~ Y_{t} = \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \beta_3 \mbox{TIME} \times \mbox{POST} + \varepsilon_t\\ \end{align} \]

\[ \begin{align} Y_t =& \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \beta_3 \mbox{POST} \times \mbox{TIME} + \\ & \beta_4 \mbox{TREAT} + \beta_5 \mbox{TREAT} \times \mbox{TIME } + \\ & \beta_6 \mbox{TREAT} \times \mbox{POST} + \beta_7 \mbox{TREAT} \times \mbox{POST} \times \mbox{TIME} + \varepsilon_t \end{align} \]

where \[ \mbox{TIME } = \mbox{ {1, ... 60}},\\ \mbox{POST } = \begin{cases} 0 \mbox{ if TIME} \le 31, \\ 1 \mbox{ if TIME > 31}\\ \end{cases},\\ \mbox{TREAT } = \begin{cases} 0 \mbox{ = Fiorina 161}, \\ 1 \mbox{ = Hadley's Hope} \end{cases} \]

\[ \begin{align} Y_t =& \beta_0 + \beta_1 \mbox{TIME} + \beta_2 \mbox{POST} + \beta_3 \mbox{POST} \times \mbox{TIME} + \\ & \beta_4 \mbox{TREAT} + \beta_5 \mbox{TREAT} \times \mbox{TIME } + \\ & \beta_6 \mbox{TREAT} \times \mbox{POST} + \beta_7 \mbox{TREAT} \times \mbox{POST} \times \mbox{TIME} + \varepsilon_t \end{align} \]

\[ \begin{align} ATT(t) = & E[Y^1 | D = 1](t) - E[Y^0 | D = 1](t) \\ E[Y^1 | D = 1](t) =& (\beta_0 + \beta_1 t + \beta_2 + \beta_3 t + \beta_4 + \beta_5 t + \beta_6 + \beta_7 t) - \\ & (\beta_0 + \beta_1 t + \beta_4 + \beta_5 t) \\ =& \beta_2 + \beta_3 t + \beta_6 + \beta_7 t \\ E[Y^0 | D = 1] =& (\beta_0 + \beta_1 t + \beta_2 + \beta_3 t) - (\beta_0 + \beta_1 t) \\ =& \beta_2 + \beta_3 t \\ E[Y^1 | D = 1] - E[Y^0 | D = 1] =& \beta_6 + \beta_7 t \end{align} \]

Model

lm_cits <- lm(y ~ t*post + t*treat + treat*post + I(treat * post * t),
             data = hhf)
## # A tibble: 8 x 5
##   term                estimate std.error statistic p.value
##   <chr>                  <dbl>     <dbl>     <dbl> <chr>  
## 1 (Intercept)          655.        9.38     69.8   <0.001 
## 2 t                      3.79      0.537     7.06  <0.001 
## 3 post                  -7.17     13.9      -0.517 0.606  
## 4 treat                -35.2      13.3      -2.65  0.009  
## 5 I(treat * post * t)   -0.530     1.13     -0.468 0.641  
## 6 t:post                -0.944     0.801    -1.18  0.241  
## 7 t:treat               -1.33      0.760    -1.75  0.082  
## 8 post:treat           -13.2      19.6      -0.673 0.502

Linear combination

summary(multcomp::glht(lm_cits, 
             linfct = matrix(c(0, 0, 0, 0, 12, 0, 0, 1), 1)))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ t * post + t * treat + treat * post + I(treat * 
##     post * t), data = hhf)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)
## 1 == 0   -19.56      22.86  -0.855    0.394
## (Adjusted p values reported -- single-step method)

Difference in differences

DiD: Intuition

ATT

\[ \begin{align} ATT =& E\left[\bar{Y^1}_{t > T_0} - \bar{Y^1}_{t \le T_0} | D = 1 \right] - E\left[\bar{Y^0}_{t > T_0} - \bar{Y^0}_{t \le T_0} | D = 1 \right] \end{align} \]

Counterfactual assumption

\[ \begin{align} E\left[\bar{Y^0}_{t > T_0} - \bar{Y^0}_{t \le T_0} | D = 1 \right] &\approx E\left[\bar{Y^0}_{t > T_0} - \bar{Y^0}_{t \le T_0} | D = 0 \right] \\ \end{align} \] Then

\[ \begin{align} ATT &\approx E\left[\bar{Y^1}_{t > T_0} - \bar{Y^1}_{t \le T_0} | D = 1 \right] - E\left[\bar{Y^0}_{t > T_0} - \bar{Y^0}_{t \le T_0} | D = 0 \right] \\ \end{align} \]

A crude DiD

Colony Pre Post Difference
Fiorina 161 (comparison) 598 691 93
Hadley’s Hope (test) 583 614 31


31-93
## [1] -62

Linear model

\[ \begin{align} Y &= \beta_0 + \beta_1 \mbox{TREAT} + \beta_2 \mbox{POST} + \\ & ~~~~~~~~~~~~~\beta_3 \mbox{TREAT} \times \mbox{POST} + \varepsilon \\ E\left[\bar{Y^1}_{t > T_0} - \bar{Y^1}_{t \le T_0} | D = 1 \right] &= [\beta_0 + \beta_1 + \beta_2 + \beta_3] - [\beta_0 + \beta_1] = \beta_2 + \beta_3 \\ E\left[\bar{Y^0}_{t > T_0} - \bar{Y^0}_{t \le T_0} | D = 0 \right] &= [\beta_0 + \beta_2] - [\beta_0] = \beta_2 \\ ATT &= [\beta_2 + \beta_3] - [\beta_2] = \beta_3 \end{align} \]

Linear model

## # A tibble: 4 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)    598.       6.44     92.9  9.90e-111
## 2 treat          -15.3      9.11     -1.67 9.67e-  2
## 3 post            92.4      9.27      9.97 2.75e- 17
## 4 treat:post     -61.1     13.1      -4.66 8.45e-  6

Linear model, econ-ish specification with unit and time fixed effects

\[ \begin{align} y_{it} &= \beta_0 + \beta_1 \mbox{TREAT} \times \mbox{POST} + \alpha_i + \gamma_t + \varepsilon_{it} \\ \end{align} \]

Linear model with unit and time fixed effects

lm_did <- lm(y ~ tx_post + t_f + colony,
                  data=hhf)
## # A tibble: 62 x 5
##    term                estimate std.error statistic  p.value
##    <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)           531.       18.8     28.3   1.23e-35
##  2 tx_postTRUE           -61.1       9.54    -6.41  2.90e- 8
##  3 t_f-29                  6.62     26.1      0.253 8.01e- 1
##  4 t_f-28                 39.6      26.1      1.52  1.35e- 1
##  5 t_f-27                 27.8      26.1      1.07  2.91e- 1
##  6 t_f-26                 23.2      26.1      0.889 3.77e- 1
##  7 t_f-25                 44.5      26.1      1.70  9.36e- 2
##  8 t_f-24                 15.7      26.1      0.602 5.50e- 1
##  9 t_f-23                 45.1      26.1      1.73  8.93e- 2
## 10 t_f-22                 52.6      26.1      2.01  4.86e- 2
## 11 t_f-21                 89.0      26.1      3.41  1.19e- 3
## 12 t_f-20                 57.0      26.1      2.18  3.31e- 2
## 13 t_f-19                 44.8      26.1      1.72  9.16e- 2
## 14 t_f-18                 84.0      26.1      3.22  2.11e- 3
## 15 t_f-17                 57.2      26.1      2.19  3.25e- 2
## 16 t_f-16                 84.6      26.1      3.24  1.97e- 3
## 17 t_f-15                 68.4      26.1      2.62  1.12e- 2
## 18 t_f-14                 95.6      26.1      3.66  5.41e- 4
## 19 t_f-13                 55.0      26.1      2.11  3.96e- 2
## 20 t_f-12                 92.7      26.1      3.55  7.66e- 4
## 21 t_f-11                 63.7      26.1      2.44  1.78e- 2
## 22 t_f-10                108.       26.1      4.12  1.22e- 4
## 23 t_f-9                  76.9      26.1      2.94  4.65e- 3
## 24 t_f-8                  60.7      26.1      2.33  2.36e- 2
## 25 t_f-7                  79.6      26.1      3.05  3.44e- 3
## 26 t_f-6                  84.5      26.1      3.24  2.00e- 3
## 27 t_f-5                  66.0      26.1      2.53  1.43e- 2
## 28 t_f-4                  91.8      26.1      3.52  8.52e- 4
## 29 t_f-3                 128.       26.1      4.92  7.63e- 6
## 30 t_f-2                 119.       26.1      4.56  2.65e- 5
## 31 t_f-1                 134.       26.1      5.12  3.63e- 6
## 32 t_f0                  102.       26.1      3.89  2.62e- 4
## 33 t_f1                  140.       26.5      5.26  2.14e- 6
## 34 t_f2                  159.       26.5      5.98  1.49e- 7
## 35 t_f3                   93.2      26.5      3.51  8.71e- 4
## 36 t_f4                  174.       26.5      6.55  1.64e- 8
## 37 t_f5                  148.       26.5      5.59  6.32e- 7
## 38 t_f6                  134.       26.5      5.05  4.63e- 6
## 39 t_f7                  142.       26.5      5.36  1.48e- 6
## 40 t_f8                  165.       26.5      6.20  6.35e- 8
## 41 t_f9                  101.       26.5      3.79  3.61e- 4
## 42 t_f10                 168.       26.5      6.31  4.12e- 8
## 43 t_f11                 135.       26.5      5.10  3.88e- 6
## 44 t_f12                 157.       26.5      5.92  1.83e- 7
## 45 t_f13                 155.       26.5      5.85  2.42e- 7
## 46 t_f14                 163.       26.5      6.12  8.49e- 8
## 47 t_f15                 178.       26.5      6.69  9.70e- 9
## 48 t_f16                 139.       26.5      5.23  2.39e- 6
## 49 t_f17                 158.       26.5      5.95  1.65e- 7
## 50 t_f18                 183.       26.5      6.90  4.31e- 9
## 51 t_f19                 210.       26.5      7.91  8.74e-11
## 52 t_f20                 147.       26.5      5.54  7.80e- 7
## 53 t_f21                 161.       26.5      6.08  1.01e- 7
## 54 t_f22                 141.       26.5      5.30  1.91e- 6
## 55 t_f23                 161.       26.5      6.08  1.01e- 7
## 56 t_f24                 194.       26.5      7.30  9.15e-10
## 57 t_f25                 177.       26.5      6.68  1.01e- 8
## 58 t_f26                 198.       26.5      7.46  4.90e-10
## 59 t_f27                 175.       26.5      6.58  1.47e- 8
## 60 t_f28                 184.       26.5      6.92  4.00e- 9
## 61 t_f29                 204.       26.5      7.67  2.20e-10
## 62 colonyHadley's Hope   -15.3       6.63    -2.30  2.50e- 2

Difference in differences

  • Assumes the pre-post difference in the control group is a reasonable surrogate for the counterfactual condition of the test group;

  • Traditionally, we would say the assumption is that the test and control groups’ longitudinal trends in the pre-intervention period are parallel. See Health Policy Data Science Lab’s website for a nuanced discussion.

  • Also assumes common shocks.

  • Accounts for baseline differences in outcomes between the two groups and shocks unrelated to the intervention that affected both groups

  • Does not account for Shocks unrelated to the intervention that only occurred in one group or differences between groups in the longitudinal trajectory before the intervention

Assumptions

  • “Parallel trends”?
  • “Counterfactual assumption”?

Parallel trends

Test for parallel trends

  • Match or weight by propensity score — but don’t include the outcome
  • Explicitly model the difference between trends — but misspecification can introduce bias
  • Choose a better control group
  • Give up. Go do a different study.

Raise your hand if you just love all of this.

Let’s conduct our analysis conditional on getting \(p > 0.05\) on a hypothesis test where the thing we’re trying to show isn’t a thing is set up as the null hypothesis.

Do you have a better idea?

  • Think more carefully about how the ATT is defined; “the test of parallel trends is neither necessary nor sufficient.”
  • Equivalence testing (Hartman and Hidalgo 2018)
  • Model specification that accounts for non-parallel trends (Bilinski and Hatfield 2019)
  • Extrapolate the deviation from the pre-intervention period into the post-intervention period (Rambachan and Roth 2019)

Fixed effects

Colonies

Three methods

  • Difference score
  • Unit and time dummy variables
  • De-meaning

Notation

\[ y_{it} = \mu_t + \beta \mathbf{x}_{it} + \gamma \mathbf{z_i} + \alpha_i + \varepsilon_{it} \] \[ \begin{align} \mu_t &~ \mbox{ is a time-varying intercept} \\ \mathbf{x}_{it} &~ \mbox{ is a vector of time-varying covariates} \\ \mathbf{z}_{i} &~ \mbox{ is a vector of time-invariant covariates} \\ \beta, \gamma &~ \mbox{ are coefficients} \\ \mathbf{\alpha}_i &~\mbox{ is time-invariant error for individuals} \\ \varepsilon_{it} &~ \mbox{ is error for each individual at each time} \end{align} \]

Two periods

\[ \begin{align} y_{i1} & = \mu_1 + \beta \mathbf{x}_{i1} + \gamma \mathbf{z}_i + \alpha_i + \varepsilon_{i1} \\ y_{i2} & = \mu_2 + \beta \mathbf{x}_{i2} + \gamma \mathbf{z}_i + \alpha_i + \varepsilon_{i2} \\ y_{i2} - y_{i1} &= (\mu_2 - \mu_1) + \beta(\mathbf{x}_{i2} - \mathbf{x}_{i1}) + (\varepsilon_{i1} - \varepsilon_{i2}) \\ \Delta y_i &= \Delta \mu + \beta \Delta \mathbf{x}_i + \Delta \varepsilon_i \end{align} \]

  • \(\alpha_i\) and \(\gamma \mathbf{z}_i\) go away
  • Cannot estimate \(\gamma\)

Two-period data, many colonies

## # A tibble: 20 x 6
##    colony        treat  post tx_post weather_index     y
##    <chr>         <dbl> <dbl>   <dbl>         <dbl> <dbl>
##  1 Arceon            1     0       0          5.36  610.
##  2 Arceon            1     1       1          5.79  644.
##  3 Arcturus          1     0       0          5.33  603.
##  4 Arcturus          1     1       1          5.14  637.
##  5 Argos             0     0       0          5.03  592.
##  6 Argos             0     1       0          4.67  679.
##  7 Articus           1     0       0          6.24  602.
##  8 Articus           1     1       1          5.02  632.
##  9 Atlas             0     0       0          5.54  562.
## 10 Atlas             0     1       0          5.41  656.
## 11 Crysalis          0     0       0          4.68  582.
## 12 Crysalis          0     1       0          5.82  683.
## 13 Cyrus             0     0       0          4.78  573.
## 14 Cyrus             0     1       0          4.80  666.
## 15 Fiorina 161       0     0       0          5.29  598.
## 16 Fiorina 161       0     1       0          4.57  691.
## 17 Hadley's Hope     1     0       0          5.88  583.
## 18 Hadley's Hope     1     1       1          4.49  614.
## 19 Relitor           1     0       0          4.52  559.
## 20 Relitor           1     1       1          4.78  607.

Calculate differences

(colonies_diff <- colonies_2p %>%
       pivot_wider(names_from = post,
                   values_from = c(tx_post, weather_index, y)) %>%
       mutate(delta_y = y_1 - y_0,
              delta_txpost = tx_post_1 - tx_post_0,
              delta_weather = weather_index_1 - weather_index_0) %>%
   select(colony, starts_with("delta")))
## # A tibble: 10 x 4
##    colony        delta_y delta_txpost delta_weather
##    <chr>           <dbl>        <dbl>         <dbl>
##  1 Arceon           34.0            1        0.437 
##  2 Arcturus         34.0            1       -0.187 
##  3 Argos            86.6            0       -0.356 
##  4 Articus          30.7            1       -1.22  
##  5 Atlas            94.3            0       -0.125 
##  6 Crysalis        101.             0        1.14  
##  7 Cyrus            93.8            0        0.0146
##  8 Fiorina 161      92.4            0       -0.714 
##  9 Hadley's Hope    31.3            1       -1.40  
## 10 Relitor          48.1            1        0.263

Difference score method

\[ \Delta y_i = \Delta \mu + \beta_1 \Delta \mbox{TREAT}_i + \beta_2 \mbox{WEATHER} + \Delta \varepsilon_i \]

lm_fe1 <- lm(delta_y ~ delta_txpost + delta_weather, data = colonies_diff)
tidy(lm_fe1)
## # A tibble: 3 x 5
##   term          estimate std.error statistic  p.value
##   <chr>            <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)      93.7       2.15     43.5  8.82e-10
## 2 delta_txpost    -55.7       3.18    -17.5  4.80e- 7
## 3 delta_weather     5.58      2.20      2.54 3.88e- 2

Unit and time dummy variables

\[ \begin{align} y_{it} & = \beta_0 + \beta_1 \mbox{TREAT} \times \mbox{POST} + \alpha_i + \gamma_t + \varepsilon_{it} , \\\\ \mbox{where} \\ \alpha_i & \mbox{ are colony fixed effects, and } \\ \gamma_t & \mbox{ are time fixed effects} \end{align} \]

Unit and time dummy variables, 2 periods

lm_fe2 <- lm(y ~ tx_post + weather_index + colony + post, data=colonies_2p)
## # A tibble: 13 x 5
##    term                estimate std.error statistic  p.value
##    <chr>                  <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)           577.       13.0      44.4  7.67e-10
##  2 tx_post               -55.7       3.18    -17.5  4.80e- 7
##  3 weather_index           5.58      2.20      2.54 3.88e- 2
##  4 colonyArcturus         -4.55      3.49     -1.31 2.33e- 1
##  5 colonyArgos           -15.1       4.25     -3.54 9.48e- 3
##  6 colonyArticus          -9.86      3.41     -2.89 2.32e- 2
##  7 colonyAtlas           -45.0       3.79    -11.9  6.85e- 6
##  8 colonyCrysalis        -20.6       3.91     -5.28 1.15e- 3
##  9 colonyCyrus           -30.7       4.32     -7.10 1.94e- 4
## 10 colonyFiorina 161      -6.37      4.17     -1.53 1.71e- 1
## 11 colonyHadley's Hope   -25.7       3.51     -7.33 1.59e- 4
## 12 colonyRelitor         -38.3       3.97     -9.65 2.71e- 5
## 13 post                   93.7       2.15     43.5  8.82e-10

Unit and time dummy variables, 60 periods

lm_fe3 <- lm(y ~ tx_post + weather_index + colony + t_f, data=colonies)
## # A tibble: 71 x 5
##    term                estimate std.error statistic   p.value
##    <chr>                  <dbl>     <dbl>     <dbl>     <dbl>
##  1 (Intercept)           530.       7.35     72.2   4.41e-276
##  2 tx_postTRUE           -56.1      3.45    -16.2   1.86e- 48
##  3 weather_index           4.80     0.315    15.2   9.52e- 44
##  4 colonyArcturus         -4.87     3.85     -1.26  2.07e-  1
##  5 colonyArgos           -15.8      4.21     -3.75  1.94e-  4
##  6 colonyArticus          -9.89     3.85     -2.57  1.05e-  2
##  7 colonyAtlas           -45.3      4.20    -10.8   1.11e- 24
##  8 colonyCrysalis        -21.2      4.20     -5.04  6.39e-  7
##  9 colonyCyrus           -31.5      4.21     -7.49  2.93e- 13
## 10 colonyFiorina 161      -7.16     4.21     -1.70  8.92e-  2
## 11 colonyHadley's Hope   -26.1      3.85     -6.78  3.26e- 11
## 12 colonyRelitor         -39.2      3.86    -10.2   2.97e- 22
## 13 t_f-29                  1.36     9.44      0.144 8.85e-  1
## 14 t_f-28                 18.3      9.45      1.93  5.38e-  2
## 15 t_f-27                 17.1      9.44      1.81  7.03e-  2
## 16 t_f-26                 11.4      9.45      1.21  2.27e-  1
## 17 t_f-25                 20.1      9.46      2.12  3.45e-  2
## 18 t_f-24                 22.0      9.44      2.33  2.00e-  2
## 19 t_f-23                 23.2      9.46      2.45  1.46e-  2
## 20 t_f-22                 30.2      9.46      3.19  1.50e-  3
## 21 t_f-21                 37.6      9.44      3.98  7.71e-  5
## 22 t_f-20                 54.0      9.44      5.72  1.76e-  8
## 23 t_f-19                 38.8      9.44      4.11  4.58e-  5
## 24 t_f-18                 44.9      9.45      4.75  2.58e-  6
## 25 t_f-17                 51.2      9.45      5.42  9.12e-  8
## 26 t_f-16                 50.0      9.44      5.30  1.74e-  7
## 27 t_f-15                 52.9      9.45      5.60  3.46e-  8
## 28 t_f-14                 58.4      9.47      6.17  1.39e-  9
## 29 t_f-13                 59.8      9.44      6.34  5.07e- 10
## 30 t_f-12                 62.7      9.44      6.64  7.70e- 11
## 31 t_f-11                 59.2      9.45      6.26  7.96e- 10
## 32 t_f-10                 69.5      9.44      7.36  7.10e- 13
## 33 t_f-9                  59.2      9.45      6.26  7.99e- 10
## 34 t_f-8                  63.6      9.46      6.72  4.75e- 11
## 35 t_f-7                  70.9      9.44      7.51  2.54e- 13
## 36 t_f-6                  82.1      9.44      8.69  4.47e- 17
## 37 t_f-5                  75.3      9.44      7.97  9.75e- 15
## 38 t_f-4                  76.8      9.44      8.14  2.85e- 15
## 39 t_f-3                  87.3      9.46      9.23  6.38e- 19
## 40 t_f-2                  99.9      9.45     10.6   8.00e- 24
## 41 t_f-1                  84.8      9.45      8.98  4.84e- 18
## 42 t_f0                   91.2      9.45      9.66  2.01e- 20
## 43 t_f1                  107.       9.61     11.1   5.07e- 26
## 44 t_f2                  119.       9.60     12.4   4.36e- 31
## 45 t_f3                  106.       9.60     11.0   1.50e- 25
## 46 t_f4                  121.       9.62     12.6   4.02e- 32
## 47 t_f5                  128.       9.60     13.3   4.41e- 35
## 48 t_f6                  122.       9.59     12.7   2.19e- 32
## 49 t_f7                  125.       9.60     13.1   4.68e- 34
## 50 t_f8                  139.       9.61     14.4   4.66e- 40
## 51 t_f9                  123.       9.59     12.8   7.68e- 33
## 52 t_f10                 145.       9.60     15.1   5.52e- 43
## 53 t_f11                 128.       9.60     13.4   2.25e- 35
## 54 t_f12                 131.       9.61     13.6   1.66e- 36
## 55 t_f13                 143.       9.60     14.9   4.13e- 42
## 56 t_f14                 131.       9.61     13.6   1.69e- 36
## 57 t_f15                 146.       9.61     15.2   1.15e- 43
## 58 t_f16                 148.       9.60     15.4   1.67e- 44
## 59 t_f17                 150.       9.61     15.6   1.40e- 45
## 60 t_f18                 152.       9.61     15.8   2.53e- 46
## 61 t_f19                 154.       9.66     15.9   6.48e- 47
## 62 t_f20                 149.       9.59     15.5   4.30e- 45
## 63 t_f21                 158.       9.61     16.4   2.38e- 49
## 64 t_f22                 143.       9.59     14.9   2.73e- 42
## 65 t_f23                 169.       9.62     17.5   1.46e- 54
## 66 t_f24                 165.       9.60     17.2   6.94e- 53
## 67 t_f25                 169.       9.61     17.6   8.10e- 55
## 68 t_f26                 185.       9.59     19.3   4.33e- 63
## 69 t_f27                 173.       9.60     18.0   8.30e- 57
## 70 t_f28                 182.       9.61     18.9   1.84e- 61
## 71 t_f29                 180.       9.61     18.7   2.66e- 60

Variance

  • Key assumption in linear regression is independent observations.
  • Repeated observations within each colony are correlated.
  • Failure to account for this correlation will usually (but not always) bias standard errors downward, yielding optimistic precision.

“STATA”-flavored standard errors

Eicker-Huber-White variance

\[ \begin{align} &\frac{N}{N-K}(\mathbf{X}^\prime \mathbf{X})^{-1} \mbox{diag} [e^2_i] \mathbf{X} (\mathbf{X}^\prime \mathbf{X})^{-1}, \\ &df = N-K \end{align} \]

In general, clustering at a higher level accounts for nested levels of correlation

vcovHC and lmtest

# Fit the model
lm_fe3 <- lm(y ~ tx_post + colony + weather_index + t_f, data=colonies)

# Estimate heteroskedasticity-consistent covariance matrix
vcv_lm_fe <- sandwich::vcovHC(lm_fe3, cluster=colonies$colony)

# Coeftest
lm_fe3_hc <- lmtest::coeftest(lm_fe3, vcv_lm_fe)

OLS versus HC standard errors

# OLS results
tidy(lm_fe3) %>% filter(!str_detect(term, "^colony|^t_f"))
## # A tibble: 3 x 5
##   term          estimate std.error statistic   p.value
##   <chr>            <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)     530.       7.35       72.2 4.41e-276
## 2 tx_postTRUE     -56.1      3.45      -16.2 1.86e- 48
## 3 weather_index     4.80     0.315      15.2 9.52e- 44
# HC results
tidy(lm_fe3_hc) %>%
  filter(!str_detect(term, "^colony|^t_f"))
## # A tibble: 3 x 5
##   term          estimate std.error statistic   p.value
##   <chr>            <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)     530.       6.26       84.7 6.75e-310
## 2 tx_postTRUE     -56.1      3.67      -15.3 6.67e- 44
## 3 weather_index     4.80     0.329      14.6 1.04e- 40

De-meaning

\[ \begin{align} y_{it} & = \beta_0 + \mathbf{\beta x}_{it} + \alpha_i + \gamma_t + \varepsilon_{it} \\ y^*_{it} & = \beta_0 + \beta \mathbf{x}^*_{it} + \gamma_t + \varepsilon_{it} , \\\\ y^*_{it} & = y_{it} - \frac{1}{n_i}\sum_t y_{it}, ~~~~~~~~~~~~~~~~\mathbf{x}^*_{it} = \mathbf{x}_{it} - \frac{1}{n_i}\sum_t \mathbf{x}_{it} \end{align} \]

Mean deviation algorithm

lm_demean <- colonies %>%
  group_by(colony) %>%
  mutate(y_star = y - mean(y),
         tx_post_star = as.numeric(tx_post) - mean(as.numeric(tx_post)),
         weather_star = weather_index - mean(weather_index)) %>%
  lm(y_star ~ tx_post_star + weather_star + t_f,
     data=.)

# Results
tidy(lmtest::coeftest(lm_demean, 
                      sandwich::vcovHC(lm_demean, cluster=colonies$colony))) %>%
  filter(!str_detect(term, "^colony|^t_f"))
## # A tibble: 3 x 5
##   term         estimate std.error statistic  p.value
##   <chr>           <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    -96.1      5.14      -18.7 2.33e-60
## 2 tx_post_star   -56.1      3.61      -15.5 3.14e-45
## 3 weather_star     4.80     0.323      14.8 6.02e-42

Package estimatr

tidy(estimatr::lm_robust(y ~ tx_post + weather_index,
          fixed_effects = ~ colony + t_f,
          clusters = colony,
          se_type = "stata",
          data=colonies))
##            term   estimate std.error statistic      p.value   conf.low
## 1   tx_postTRUE -56.066789 2.8722326 -19.52028 1.124640e-08 -62.564231
## 2 weather_index   4.795408 0.2413074  19.87261 9.605851e-09   4.249533
##    conf.high df outcome
## 1 -49.569348  9       y
## 2   5.341283  9       y

Benchmarking

bench <- map_dfr(1:100, 
                 function(i) {
                   tic(paste0("lm() iteration ", i))
                   # Fit the model
                   lm_fe <- lm(y ~ tx_post + weather_index + colony + t_f,
                               data=colonies)
                   
                   # Estimate heteroskedasticity-consistent covariance matrix
                   vcv_lm_fe <- sandwich::vcovHC(lm_fe, cluster=colonies$colony) 
                   
                   # Results
                   rslts <- tidy(lmtest::coeftest(lm_fe, vcv_lm_fe)) %>%
                     filter(!str_detect(term, "^colony|^t_f"))
                   tt_lm <- toc()
                   
                   tic(paste0("lm_robust() iteration ", i))
                   rslts2 <- tidy(estimatr::lm_robust(y ~ tx_post + weather_index,
                                                      fixed_effects = ~ colony + t_f,
                                                      clusters = colony,
                                                      se_type = "stata",
                                                      data=colonies))
                   tt_lmr <- toc()
                   
                   return(tibble(lm = tt_lm$toc - tt_lm$tic,
                                 lm_robust = tt_lmr$toc - tt_lmr$tic))
                 })
## lm() iteration 1: 0.035 sec elapsed
## lm_robust() iteration 1: 0.301 sec elapsed
## lm() iteration 2: 0.019 sec elapsed
## lm_robust() iteration 2: 0.008 sec elapsed
## lm() iteration 3: 0.03 sec elapsed
## lm_robust() iteration 3: 0.01 sec elapsed
## lm() iteration 4: 0.021 sec elapsed
## lm_robust() iteration 4: 0.009 sec elapsed
## lm() iteration 5: 0.033 sec elapsed
## lm_robust() iteration 5: 0.009 sec elapsed
## lm() iteration 6: 0.022 sec elapsed
## lm_robust() iteration 6: 0.008 sec elapsed
## lm() iteration 7: 0.034 sec elapsed
## lm_robust() iteration 7: 0.008 sec elapsed
## lm() iteration 8: 0.031 sec elapsed
## lm_robust() iteration 8: 0.008 sec elapsed
## lm() iteration 9: 0.033 sec elapsed
## lm_robust() iteration 9: 0.009 sec elapsed
## lm() iteration 10: 0.022 sec elapsed
## lm_robust() iteration 10: 0.008 sec elapsed
## lm() iteration 11: 0.033 sec elapsed
## lm_robust() iteration 11: 0.008 sec elapsed
## lm() iteration 12: 0.022 sec elapsed
## lm_robust() iteration 12: 0.008 sec elapsed
## lm() iteration 13: 0.037 sec elapsed
## lm_robust() iteration 13: 0.009 sec elapsed
## lm() iteration 14: 0.021 sec elapsed
## lm_robust() iteration 14: 0.029 sec elapsed
## lm() iteration 15: 0.019 sec elapsed
## lm_robust() iteration 15: 0.008 sec elapsed
## lm() iteration 16: 0.021 sec elapsed
## lm_robust() iteration 16: 0.01 sec elapsed
## lm() iteration 17: 0.031 sec elapsed
## lm_robust() iteration 17: 0.008 sec elapsed
## lm() iteration 18: 0.021 sec elapsed
## lm_robust() iteration 18: 0.007 sec elapsed
## lm() iteration 19: 0.031 sec elapsed
## lm_robust() iteration 19: 0.008 sec elapsed
## lm() iteration 20: 0.021 sec elapsed
## lm_robust() iteration 20: 0.008 sec elapsed
## lm() iteration 21: 0.041 sec elapsed
## lm_robust() iteration 21: 0.007 sec elapsed
## lm() iteration 22: 0.02 sec elapsed
## lm_robust() iteration 22: 0.008 sec elapsed
## lm() iteration 23: 0.031 sec elapsed
## lm_robust() iteration 23: 0.008 sec elapsed
## lm() iteration 24: 0.021 sec elapsed
## lm_robust() iteration 24: 0.008 sec elapsed
## lm() iteration 25: 0.033 sec elapsed
## lm_robust() iteration 25: 0.008 sec elapsed
## lm() iteration 26: 0.022 sec elapsed
## lm_robust() iteration 26: 0.008 sec elapsed
## lm() iteration 27: 0.042 sec elapsed
## lm_robust() iteration 27: 0.007 sec elapsed
## lm() iteration 28: 0.021 sec elapsed
## lm_robust() iteration 28: 0.008 sec elapsed
## lm() iteration 29: 0.032 sec elapsed
## lm_robust() iteration 29: 0.008 sec elapsed
## lm() iteration 30: 0.02 sec elapsed
## lm_robust() iteration 30: 0.008 sec elapsed
## lm() iteration 31: 0.03 sec elapsed
## lm_robust() iteration 31: 0.011 sec elapsed
## lm() iteration 32: 0.02 sec elapsed
## lm_robust() iteration 32: 0.009 sec elapsed
## lm() iteration 33: 0.043 sec elapsed
## lm_robust() iteration 33: 0.009 sec elapsed
## lm() iteration 34: 0.021 sec elapsed
## lm_robust() iteration 34: 0.009 sec elapsed
## lm() iteration 35: 0.028 sec elapsed
## lm_robust() iteration 35: 0.011 sec elapsed
## lm() iteration 36: 0.02 sec elapsed
## lm_robust() iteration 36: 0.008 sec elapsed
## lm() iteration 37: 0.033 sec elapsed
## lm_robust() iteration 37: 0.007 sec elapsed
## lm() iteration 38: 0.022 sec elapsed
## lm_robust() iteration 38: 0.008 sec elapsed
## lm() iteration 39: 0.041 sec elapsed
## lm_robust() iteration 39: 0.007 sec elapsed
## lm() iteration 40: 0.021 sec elapsed
## lm_robust() iteration 40: 0.008 sec elapsed
## lm() iteration 41: 0.033 sec elapsed
## lm_robust() iteration 41: 0.009 sec elapsed
## lm() iteration 42: 0.021 sec elapsed
## lm_robust() iteration 42: 0.008 sec elapsed
## lm() iteration 43: 0.029 sec elapsed
## lm_robust() iteration 43: 0.011 sec elapsed
## lm() iteration 44: 0.021 sec elapsed
## lm_robust() iteration 44: 0.02 sec elapsed
## lm() iteration 45: 0.03 sec elapsed
## lm_robust() iteration 45: 0.009 sec elapsed
## lm() iteration 46: 0.021 sec elapsed
## lm_robust() iteration 46: 0.008 sec elapsed
## lm() iteration 47: 0.032 sec elapsed
## lm_robust() iteration 47: 0.008 sec elapsed
## lm() iteration 48: 0.021 sec elapsed
## lm_robust() iteration 48: 0.008 sec elapsed
## lm() iteration 49: 0.033 sec elapsed
## lm_robust() iteration 49: 0.008 sec elapsed
## lm() iteration 50: 0.023 sec elapsed
## lm_robust() iteration 50: 0.008 sec elapsed
## lm() iteration 51: 0.042 sec elapsed
## lm_robust() iteration 51: 0.008 sec elapsed
## lm() iteration 52: 0.02 sec elapsed
## lm_robust() iteration 52: 0.008 sec elapsed
## lm() iteration 53: 0.031 sec elapsed
## lm_robust() iteration 53: 0.008 sec elapsed
## lm() iteration 54: 0.021 sec elapsed
## lm_robust() iteration 54: 0.007 sec elapsed
## lm() iteration 55: 0.032 sec elapsed
## lm_robust() iteration 55: 0.009 sec elapsed
## lm() iteration 56: 0.021 sec elapsed
## lm_robust() iteration 56: 0.008 sec elapsed
## lm() iteration 57: 0.027 sec elapsed
## lm_robust() iteration 57: 0.026 sec elapsed
## lm() iteration 58: 0.019 sec elapsed
## lm_robust() iteration 58: 0.008 sec elapsed
## lm() iteration 59: 0.03 sec elapsed
## lm_robust() iteration 59: 0.011 sec elapsed
## lm() iteration 60: 0.023 sec elapsed
## lm_robust() iteration 60: 0.008 sec elapsed
## lm() iteration 61: 0.033 sec elapsed
## lm_robust() iteration 61: 0.008 sec elapsed
## lm() iteration 62: 0.022 sec elapsed
## lm_robust() iteration 62: 0.008 sec elapsed
## lm() iteration 63: 0.023 sec elapsed
## lm_robust() iteration 63: 0.017 sec elapsed
## lm() iteration 64: 0.03 sec elapsed
## lm_robust() iteration 64: 0.008 sec elapsed
## lm() iteration 65: 0.033 sec elapsed
## lm_robust() iteration 65: 0.009 sec elapsed
## lm() iteration 66: 0.021 sec elapsed
## lm_robust() iteration 66: 0.009 sec elapsed
## lm() iteration 67: 0.034 sec elapsed
## lm_robust() iteration 67: 0.009 sec elapsed
## lm() iteration 68: 0.02 sec elapsed
## lm_robust() iteration 68: 0.008 sec elapsed
## lm() iteration 69: 0.032 sec elapsed
## lm_robust() iteration 69: 0.008 sec elapsed
## lm() iteration 70: 0.03 sec elapsed
## lm_robust() iteration 70: 0.008 sec elapsed
## lm() iteration 71: 0.031 sec elapsed
## lm_robust() iteration 71: 0.008 sec elapsed
## lm() iteration 72: 0.02 sec elapsed
## lm_robust() iteration 72: 0.008 sec elapsed
## lm() iteration 73: 0.032 sec elapsed
## lm_robust() iteration 73: 0.009 sec elapsed
## lm() iteration 74: 0.022 sec elapsed
## lm_robust() iteration 74: 0.008 sec elapsed
## lm() iteration 75: 0.031 sec elapsed
## lm_robust() iteration 75: 0.011 sec elapsed
## lm() iteration 76: 0.029 sec elapsed
## lm_robust() iteration 76: 0.008 sec elapsed
## lm() iteration 77: 0.031 sec elapsed
## lm_robust() iteration 77: 0.009 sec elapsed
## lm() iteration 78: 0.021 sec elapsed
## lm_robust() iteration 78: 0.008 sec elapsed
## lm() iteration 79: 0.031 sec elapsed
## lm_robust() iteration 79: 0.009 sec elapsed
## lm() iteration 80: 0.02 sec elapsed
## lm_robust() iteration 80: 0.008 sec elapsed
## lm() iteration 81: 0.021 sec elapsed
## lm_robust() iteration 81: 0.02 sec elapsed
## lm() iteration 82: 0.031 sec elapsed
## lm_robust() iteration 82: 0.008 sec elapsed
## lm() iteration 83: 0.029 sec elapsed
## lm_robust() iteration 83: 0.011 sec elapsed
## lm() iteration 84: 0.02 sec elapsed
## lm_robust() iteration 84: 0.008 sec elapsed
## lm() iteration 85: 0.032 sec elapsed
## lm_robust() iteration 85: 0.008 sec elapsed
## lm() iteration 86: 0.022 sec elapsed
## lm_robust() iteration 86: 0.009 sec elapsed
## lm() iteration 87: 0.045 sec elapsed
## lm_robust() iteration 87: 0.008 sec elapsed
## lm() iteration 88: 0.02 sec elapsed
## lm_robust() iteration 88: 0.009 sec elapsed
## lm() iteration 89: 0.033 sec elapsed
## lm_robust() iteration 89: 0.009 sec elapsed
## lm() iteration 90: 0.02 sec elapsed
## lm_robust() iteration 90: 0.008 sec elapsed
## lm() iteration 91: 0.03 sec elapsed
## lm_robust() iteration 91: 0.01 sec elapsed
## lm() iteration 92: 0.021 sec elapsed
## lm_robust() iteration 92: 0.008 sec elapsed
## lm() iteration 93: 0.043 sec elapsed
## lm_robust() iteration 93: 0.009 sec elapsed
## lm() iteration 94: 0.021 sec elapsed
## lm_robust() iteration 94: 0.008 sec elapsed
## lm() iteration 95: 0.032 sec elapsed
## lm_robust() iteration 95: 0.008 sec elapsed
## lm() iteration 96: 0.021 sec elapsed
## lm_robust() iteration 96: 0.008 sec elapsed
## lm() iteration 97: 0.032 sec elapsed
## lm_robust() iteration 97: 0.009 sec elapsed
## lm() iteration 98: 0.032 sec elapsed
## lm_robust() iteration 98: 0.009 sec elapsed
## lm() iteration 99: 0.031 sec elapsed
## lm_robust() iteration 99: 0.008 sec elapsed
## lm() iteration 100: 0.021 sec elapsed
## lm_robust() iteration 100: 0.008 sec elapsed

Benchmarking

with(bench, t.test(lm, lm_robust))
## 
##  Welch Two Sample t-test
## 
## data:  lm and lm_robust
## t = 5.0652, df = 109.49, p-value = 1.665e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.009295288 0.021244712
## sample estimates:
## mean of x mean of y 
##   0.02732   0.01205

Comparison to mixed-effects models

Key difference:

  • Mixed effects models assume random effects are uncorrelated covariates;
  • Fixed effects models assume fixed effects are perfectly correlated with time-invariant covariates, both measured and unmeasured.
    • Cannot estimate coefficients for time-invariant covariates
    • Remove potential bias from unmeasured time-invariant confounders

Comparison to mixed-effects models

  • Fixed effects are not constrained in distribution, whereas random effects are \(\sim N(0, \sigma_b)\).
  • Fixed-effects models base estimates of \({\bf \hat{\beta}}\) exclusively on within-subject variance, adjusting out all of the between-subject variance, whereas mixed-effects models utilize both sources of variance, which can yield much smaller standard errors.
  • Bias-variance trade-off.

Linear probability models

Controversy

“There’s nothing wrong with OLS.” — John McConnell

“OLS is frequently a biased estimator and almost always an inconsistent estimator of the LPM.” — Horrace and Oaxaca

“The LPM won’t give the true marginal effects from the right nonlinear model. But then, the same is true for the “wrong” nonlinear model! The fact that we have a probit, a logit, and the LPM is just a statement to the fact that we don’t know what the “right” model is. Hence, there is a lot to be said for sticking to a linear regression function as compared to a fairly arbitrary choice of a non-linear one!" — Steffen Pischke

Practicalities

  • A link function can dramatically increase the complexity of a difference-in-differences model.
  • In my own experiments, predicted probability from a logistic model has always been within a tiny margin of the prediction from the corresponding LPM.
    • But I always check.
  • Bias increases as the proportion of predicted probabilities outside [0,1] increases.
  • Consider the PI, the target journal, and the audience.

Risk of alien infestation

LPM

tidy(lpm_infest <- lm(infest ~ treat*post + t + weather_index, 
                      data=aliens))
## # A tibble: 6 x 5
##   term          estimate std.error statistic  p.value
##   <chr>            <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)     0.248    0.0518      4.80  2.02e- 6
## 2 treat           0.0166   0.0449      0.370 7.11e- 1
## 3 post            0.205    0.0721      2.85  4.52e- 3
## 4 t               0.0131   0.00186     7.03  5.59e-12
## 5 weather_index   0.0228   0.00562     4.07  5.44e- 5
## 6 treat:post     -0.465    0.0645     -7.20  1.79e-12

Logit

tidy(glm_infest <- glm(infest ~ treat*post + t + weather_index, 
                       data=aliens,
                       family=binomial(link=logit)),
     exponentiate=TRUE)
## # A tibble: 6 x 5
##   term          estimate std.error statistic  p.value
##   <chr>            <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)     0.272     0.339     -3.84  1.21e- 4
## 2 treat           1.10      0.313      0.293 7.69e- 1
## 3 post            1.86      0.431      1.44  1.50e- 1
## 4 t               1.09      0.0131     6.54  6.02e-11
## 5 weather_index   1.16      0.0380     3.96  7.41e- 5
## 6 treat:post      0.0946    0.427     -5.53  3.27e- 8

Predicted probabilities

Fun stuff to read

  1. Health Policy Data Science Lab’s Difference-in-Differences page
  2. Lindner S, McConnell KJ. Difference-in-differences and matching on outcomes: a tale of two unobservables. Health Serv Outcomes Res Method. 2019;19(2-3):127-144. doi:10.1007/s10742-018-0189-0
  3. Bilinski A, Hatfield LA. Nothing to see here? Non-inferiority approaches to parallel trends and other model assumptions. arXiv:180503273 [stat]. Published online January 16, 2020. Accessed April 30, 2020. http://arxiv.org/abs/1805.03273
  4. Rambachan A, Roth J. An Honest Approach to Parallel Trends. Working paper.
  5. Revisiting the Difference-in-Differences Parallel Trends Assumption: Part I Pre-Trend Testing. Accessed November 13, 2020. https://blogs.worldbank.org/impactevaluations/revisiting-difference-differences-parallel-trends-assumption-part-i-pre-trend

More fun stuff to read

  1. Revisiting the Difference-in-Differences Parallel Trends Assumption: Part II What happens if the parallel trends assumption is (might be) violated? Accessed November 13, 2020. https://blogs.worldbank.org/impactevaluations/revisiting-difference-differences-parallel-trends-assumption-part-ii-what-happens
  2. Horrace WC, Oaxaca RL. Results on the bias and inconsistency of ordinary least squares for the linear probability model. Economics Letters. 2006;90(3):321-327. doi:10.1016/j.econlet.2005.08.024
  3. Whether to probit or to probe it: in defense of the Linear Probability Model. Accessed November 13, 2020. https://blogs.worldbank.org/impactevaluations/whether-to-probit-or-to-probe-it-in-defense-of-the-linear-probability-model

.